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Motivated by recent experiments on generation of wave patterns by a polariton condensate incident 
on a localized obstacle, we study the characteristics of such flows under the condition that irreversible 
processes play a crucial role in the system. The dynamics of a non-resonantly pumped polariton 
condensate in a quasi-one-dimensional quantum wire is modeled by a Gross-Pitaevskii equation 
with additional phenomenological terms accounting for the dissipation and pumping processes. The 
response of the condensate flow to an external potential describing a localized obstacle is considered 
in the weak-perturbation limit and also in the nonlinear regime. The transition from a viscous drag 
to a regime of wave resistance is identifled and studied in detail. 

PACS numbers: 03.75.Kk,71.36.+c 



C/3 

bJQf 



o 



oo 
o 



X: 



I. INTRODUCTION 

The ability to move with respect to an obstacle with- 
out dissipating energy is one of the most intuitive and ap- 
pealing definition of superfluidity. This is the reason why 
the motion of quantum fluids with respect to obstacles 
has been used in several experiments aiming at revealing 
a superfluid behavior in different physical systems: *He 
(see, e.g., Refs. [l| and|3), '^He (Ref. ultracold atomic 
vapors^""— and more recently polariton condensates^rJ^. 

For a weakly perturbing impurity moving at constant 
velocity V in a conservative atomic Bose-Einstein con- 
densed (BEC) system at zero temperature, the Lan- 
dau criterioni^ predicts that there exists a critical ve- 
locity 14rit separating two different behaviors: (i) for 
V < Vcvit no excitations are emitted away from the ob- 
stacle and, hence, there is no drag force; (ii) for V > Vait 
a Cherenkov radiation of linear waves occurs; these waves 
carry momentum away from the impurity which is thus 
subject to a finite drag force. The first regime is super- 
fiuid and the second one is dissipativei^. 

In a pumped non-equilibrium polariton condensate, 
even when kinematically allowed, propagating distur- 
bances are always damped due to the finite lifetime of 
the polaritons. As a result, the well defined transition 
between superfiuid and dissipative regimes transforms 
in these damped systems into a crossover characterized 
by different forms of wave patterns: localized for small 
enough flow velocity; oscillatory and extended for large 
enough flow velocity. The boundary between these two 
regimes is typically not abrupt: just at the transition 
point the decay length of a propagating wave is less than 
its wavelength and this disturbance can hardly be distin- 
guished from a localized perturbation. It might thus be 
difficult to separate a superfluid regime from a dissipative 
one by studying the wave pattern created by an obstacle. 
Nevertheless, the concept of superfluidity is often em- 
ployed because it permits a simple qualitative discussion 
of the processes taking place in the flow of a polariton 



condensate. 

In the present work we study in detail the wake of a 
polariton condensate past an obstacle and the associated 
drag force. We argue that, for low enough damping, the 
superfluid/dissipative transition is better understood in 
term of a crossover of the force experienced by the ob- 
stacle from a viscous drag to wave resistance, in analogy 
to what is observed for capillarity-gravity waves. 

The paper is organized as follows. In Sec.|lT]we present 
the phenomenological one-dimensional model we use and 
present our strategy for studying the speciflc features of 
typical flows. In Sec. IIIII we set up a general pertur- 
bative analysis of the motion of the polariton gas past 
a weak obstacle and discuss the domain of validity of 
this approach. In Sec. II VI we obtain non-perturbative re- 
sults valid for a localized narrow impurity using several 
approximation schemes (the so-called hydraulic approxi- 
mation in Sec. IIV Al and Whitham averaging method in 
Sec. IIVB[) and also numerical integration (Sec. IIV CI) . 
Finally we present our conclusions in Sec. [V] Some tech- 
nical points are given in the appendices. In Appendix \X\ 
we study the poles of the response function of the system 
and in Appendix [B] we present the Whitham theory we 
use in Sec. lIVBl of the main text. 



II. THE MODEL 

We study the flow of a polariton condensate past an ob- 
stacle disregarding possible effects of polarization of the 
light modes in the cavity. We consider a configuration in 
which excitons are confined in a one-dimensional quan- 
tum wire and, as a result, the polariton condensate is 
described by an order parameter ^p{x, t) whose dynamics 
is modeled by a Gross-Pitaevskii equation of the form 

ihipt = -—il^xx + {Uc^t{x,t) + ap)^p + i{'y-Tp)^. (1) 

In Eq. ^ m is the polariton effective mass (in the 
parabolic dispersion approximation, valid at small mo- 
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menta), p{x,t) = \ilj{x,t)\'^ is the polariton density and 
Ucxt{x,t) describes the potential of a locahzed obstacle, 
possibly in motion relative to the polariton gas. Inter- 
action effects are described by an effective local repul- 
sive term characterized by the nonlinear coupling con- 
stant a > 0. There is a whole body of evidence showing 
that the overall effective interaction between polaritons 
is repulsive. Some of the most direct manifestations of 
this repulsion are the observed emission blueshifti^"— 
and the expulsion of the condensate from a pumping 
regioi*i2i2S. Another consequence of repulsion, very im- 
portant for the present study, is the absence of scattering 
from a defect — first observed in Refs. and[l3 — and the 
related emission of nonlinear excitations^ (solitons and 
vortices) whose generation is typically associated to a loss 
of superfluidityii"— 

Due to the finite lifetime of the polaritons, the system 
needs to be pumped. Following Refs. [23 - |25l . we schemat- 
ically describe this effect by the last term of Eq. ([T|): the 
term hipt = 7 "0 phenomenologically describes the com- 
bined effects of the pumping and decay processes and 
for 7 > an overall gain leads, if not compensated, to 
an exponential increase of the density. This increase is 
counterbalanced by the term Hipt = —Vpij) (where F > 0) 
which accounts for a saturation of the gain at large den- 
sity and allows to reach a steady state configuration — 
resulting from dynamical equilibrium between gain and 
losses — with a finite density po = 7/1"- Eq. ([1]) corre- 
sponds to a situation where the pumping extends over 
all space. This models a system where an obstacle is 
present within a large reservoir, and simplifies the the- 
oretical treatment because the stationary density in ab- 
sence of external potential is constant. Results where the 
obstacle is present outside of the pumping region will be 
presented in a forthcoming publication^^. 

Localized structural defects are naturally present in 
many samples; they can also be artificially created by 
means of lithographic techniques or by a continuous- wave 
laser. If an obstacle is introduced into the condensate, the 
state with uniform density pq is disturbed. We suppose 
that the obstacle is described by a potential Uc-x.t{x,t) 
with a finite spatial extension [verifying Uc^t{x,t) — >■ 
as |a;| — >■ oo]. In many experiments the condensate is 
put into motion with respect to the obstacle by resonant 
pumping. Here we rather describe a situation with non- 
resonant pumping, where condensation can be forced to 
occur in a finite-momentum state by seeding the system 
with a short coherent-light pulse^^. However, we believe 
that the gross features of the theoretical study of the wave 
patterns and of the drag force are not essentially affected 
by the technique used for setting the fluid into motion. 
This is supported by a comparison of the results of the 
present work with the one of Ref. [13 where a continuous 
transition at a critical velocity (possibly different from 
the speed of sound) is also observed in a perturbative 
study of a resonantly driven polariton fluid. 

As just discussed, in typical experiments with polari- 
ton condensates the obstacle does not move and instead 



the condensate is put into motion with some velocity V. 
However we shall sometimes use for convenience a refer- 
ence frame in which the condensate is at rest (far enough 
from the obstacle) and where the obstacle moves with 
velocity —V: Ucy^t{x,t) = /oxt(2^ + Vt). A comprehen- 
sive study of this problem can be done in the case of an 
obstacle represented by a weak potential which induces 
a wave disturbance corresponding to small modifications 
of the parameters of the flow. In this conflguration the 
problem can be treated in the framework of perturbation 
theory which is presented in the next section of the paper. 
This corresponds to the extension to damped systems of 
previous perturbative studies of BEC atomic vapors^Sr^. 
In this approach the case of a 6-peak impurity is of spe- 
cial interest because its solution corresponds to the Green 
function of the problem, and we treat it with special care. 
In Sec. IIVI instead, we consider the wave pattern gener- 
ated by the flow of a polariton condensate past a strong 
obstacle potential, when perturbation theory does not 
longer apply. In this case, it is appropriate to distin- 
guish between wide and narrow obstacles depending on 
the ratio of their sizes to the healing length ^ (^ is the de 
Broglie wavelength of polaritons moving with the sound 
velocity; see its definition in the next paragraph). When 
a narrow obstacle moves at supersonic speed the down- 
stream profile has a rather smooth behavior which can 
be described by a dispersionless approach, the hydraulic 
approximation which we present in Sec. IIV Al On the 
other hand, the upstream-wave structure can be repre- 
sented (for small enough damping coefficient) as a weakly 
modulated nonlinear periodic wave which is a damped 
dispersive shock wave. Such shocks have been studied 
for the case of a wide obstacle with the use of Whitham 
modulation theory in Ref. |32| . In the present work we 
present a similar and more detailed study in the case of 
a (5-impurity in Sec. IIVBI 

In absence of external potential, a homogeneous and 
stationary solution of Eq. ^ corresponds to an order 
parameter of the form il]{x,t) = ^/po exp{—ipt / h) , where 
Pq is the uniform density and fi is the chemical potential. 
Inserting this expression in ([T} one finds po — ^ /V (nec- 
essary for obtaining a real p corresponding to a time- 
independent density) and p = apQ. The characteristic 
density po and energy p are associated to characteris- 
tic vel ocity an d distance, namely the speed of sound^ 
Cs = \J apo I m and the healing length £, = h/ (mcs). 

We will see below that, for a given obstacle potential 
Uc:xtix,t), the flow pattern is monitored by only two di- 
mensionless parameters: the Mach number M and the 
damping parameter rj defined as 

V 7 
M = — and 77=-. (2) 

Cs p 

Having identified the relevant parameters of the problem 
one can simplify the notations by expressing densities in 
units of Pq, distances in units of ^, times in units of ^/cg 
and energies in units of p. In these new variables Eq. ([T]) 
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takes the form 



In the case where rj is not zero one geti 



.22 



iipt = -^V'xx + {Ucxt{x,t) + p)ip + i7/(l - p)tjj. (3) 

From now on, we shah use this dimensionless form of the 
damped Gross-Pitaevskii equation. 



III. FLOW PAST A WEAK OBSTACLE 
A. General linear theory 

In absence of external potential Eq. ([3]) admits a uni- 
form stationary solution of the form ip{x,t) — exp(— it). 
If the potential of the obstacle is weak, one can evaluate 
the density and the flow velocity profiles of the polariton 
condensate perturbatively. In this case one looks for a 
solution of Eq. ([3]) of the form 



i){x,t) = [1 + ^{x^t)] exp(-it). 



(4) 



assuming that |</3(a::,t)| <C 1. Linearizing Eq. ([3]) with res- 
pect to </?(a;, t) and Uc^^ti^, t) and introducing the Fourier 
transforms 



viq,uj) 



dx dt 



R2 



ip*{x,t) 



-i{qx—ujt) 



(5) 



one finds that u{q,uj) and v{q^uj) satisfy the following 
linear system: 



where 



£ ^ ih — uj + \ - \ri 1 -irj 



1 + i?/ \+uj + l + h-jj 



(6) 



(7) 



When t/oxt ('?,<^) = 0, i.e., in the absence of the external 
obstacle, non-trivial solutions and v{q,uj) of the 

2x2 system © exist only when the determinant 



D{q,io)^q^ 



21770; 



(8) 



of the matrix £ is identically null. The resolution of the 
characteristic equation D{q,uj) = yields the dispersion 
relation uj{q) of the elementary excitations propagating 
on top of a homogeneous and stationary profile. Let us 
first consider the case 77 — >■ (and also in dimensional 
units r ^ in such a way that the density po = j/T is 
kept constant). In this case one finds that the excitation 
spectrum is the Bogoliubov one, i.e., one recovers the 
dispersion relation of elementary excitations of a weakly 
interacting atomic Bose gas: uj{q) = ±LUB{q), where 



WB(g) 



1 + 



(9) 



uj{q) = 



-i?/± 1^772 -wKq) if \q\<q* 



-b] ± V^K?) - 77^ if \q\ > q* 



where 



q,= 2 + 



1/2 



(10) 



(11) 



In the ideal case (t? = and then = 0) long- wavelength 
perturbations {\q\ <C 1) correspond to sound waves with 
a linear dispersion ujb{q) — Q and with a sound veloc- 
ity equal to unity in our dimensionless units. As an- 
nounced in note 33, perturbations with \q\ < q* do not 
propagate in presence of finite damping (77 7^ 0). How- 
ever, for small 77 there exists a finite region of wavenum- 
ber (g* <^ \q\ <C 1) for which the dispersion relation 
(jlOP can be approximated by the long-wavelength limit 
Lo{q) ~ q — 177 describing weakly damped sound-waves. 

Let us now consider the general case where Uc^t{x,t) 
is not zero: the linear waves are generated by the exter- 
nal potential and their Fourier components u{q, w) and 
v{q,uj) can be expressed by means of Eq. ^ in terms 
of this potential. This yields the following expression for 
the first order density modulation 6p = ip + ip* induced 
by C/oxt(a;,i): 



Sp{x,t) 



dq duj 



x{q,u;)UeAq,^)e'^'^-^'\ (12) 



where 



Spill ^) 
C/cxt(9,w) 



D{q,u) 



(13) 



is the linear response function of the system. A con- 
figuration of great experimental interest corresponds to 
the case where the condensate moves at constant veloc- 
ity with respect to a static obstacle. In this case, in the 
frame where the condensate is at rest, the external po- 
tential is of the form 



(14) 



where M is, in our dimensionless units, the velocity of 
the obstacle with respect to the condensate. For being 
specific, we shall henceforth consider the case M > 
which corresponds to an obstacle moving to the left 
in a frame where the condensate is at rest. Denot- 
ing by /cxt the Fourier transform of /ext [i-e., /ext(g) = 
/jj dz /ext(^) exp(— igz)] the expression of 5p{x,t) in the 
case of an external potential of the form ((T4)) reads 



5p{x, t) 



dq 
27r 



x(<z,-Mq)/ext(g)e'«(-+^'^*) 



(15) 



= / dzK{x + Mt- z)U^t{z), 
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where 



K{X) = 



ZTT 



(16) 



One can first remark that 5p is a function oi x + Mt 
only: the perturbative approach predicts that the den- 
sity modulations induced by an obstacle moving at con- 
stant velocity are stationary in the reference frame where 
the obstacle is at rest. Note however that, in absence of 
damping, experiment performed on atomic condensates^ 
and theory22ii2£!^ show that there is a regime of time- 
dependent flows for impurity velocities close to the speed 
of sound. This is a nonlinear effect which is missed by 
the perturbative approach. In presence of damping this 
time-dependent behavior also exits but, in a numerical 
study of nonlinear effects in presence of a wide obsta- 
cle, it is observed in a smaller domain in the parameter 
space (Intensity of Ue^t,V) than when r] = 0'^^. This is 
confirmed in the case of a narrow obstacle by the numer- 
ical results of Sec. I IV CI below. In this respect, the per- 
turbative result — being stationary — is thus more sound 
in presence of damping since in this case the domain of 
time-dependent flows is reduced. We make this discus- 
sion quantitative at the end of Sec. IIIIBI by discussing 
the parameters governing the mathematical validity of 
perturbation theory. 

A particular property of solution (|15p comes from the 
conservation equation 



Pt+jx = 2r;/?(l - p), 



(17) 



where j — lnl{^p*'lpx) is the particle current-density. Ac- 
tually, Eq. (IT7)) is a bona fide conservation equation only 
when 77 = 0. For nonzero 77, the number of particles 
is not conserved and Eq. (|17p should rather be called a 
"non-conservation" equation for the current of particles. 
Eq. PT)) is a direct consequence of (H}; in the stationary 
regime it implies^^i^ 



/ dxp(l -p) = 0. 

JR 



(18) 



At the perturbative level this reads dx 6p — which 
is trivially verified by (jlSp . 



B. Flow past a 5-impurity 

It is instructive to discuss in greater details the char- 
acteristics of the wave pattern induced by a localized ob- 
stacle with the potential 



Ucxtix,t) ^ >cd{x + Mt). 



Then one gets 



Sp{x,t) = m:K{X ^x + Mt). 



(19) 



(20) 



This density modulation is typical for the perturbations 
induced by a narrow obstacle moving in the polariton 
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FIG. 1. (Color online) Location of the three poles qi, q2 
and ga of xil, —Mq) in the complex q-plane. For positive 
(negative) X = x + Mt the integral in Eq. (|16|) is evaluated by 
closing the contour from above (below). As a result, for M > 
Merit (damped) density oscillations are observed upstream the 
obstacle (i.e., for X < 0). 



condensate. Besides, the solution of the (5-impurity prob- 
lem is particularly interesting because K{X) is the Green 
function from which the result for any potential is ob- 
tained by convolution [cf. Eq. (flS)) ]. 

The integral ([T6| can be computed by the method of 
residues and K{X) has different behaviors depending on 
the value of M and corresponding to different arrange- 
ments of the poles of x{Qj —Mq) in the complex g-plane. 
The poles are the roots of the equation D{q, ~Mq) /q — 
which reads 



q 



4(1 - M^)q + Bi-qM = 0. 



(21) 



The explicit expression of the three poles qi , q2 and q^ in 
function of 77 and M is given in Appendix [X] One obtains 
the following generic expression 

3 

K{X) - i ^ sgn(Im q,) Res{qi) e[sgn(Im qt)X] e'«^^, 
1=1 

(22) 

where is the Heaviside step function and Res((7^) is the 
residue of x(q, —Mq) at qt. 



Res(g£) 



-Aqi 



4(1 - M2) 



(23) 



There exists a critical velocity Merit below which the 
poles of xili —Mq) are all located on the imaginary axis 
(cf. Fig. [1] and also Appendix |^ and in this case for- 
mula shows that K{X) exponentially goes to when 
|X| — > 00. A more transparent expression can be ob- 
tained by explicitly solving the third order equation (|2ip . 
This yields 



K{X sC 0) 



K{X > 0) 



A-B 
A-3b' 



2 A + B 
A A + W' 



■AA-B)X 



AAB 

' A^ - 

^{A+B)X 



„2BX 



(24) 
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where A and B are positive real numbers {A > B ^ 0) 
depending on M and 77, whose expressions are given in 
Appendix El [Eq. (|A3|)]. 

On the other hand, when M > Merit, two of the poles 
acquire a real part and are symmetrically disposed with 
respect to the imaginary axis (cf. Fig.[T]). In this case the 
wave pattern is given by the explicit formulas 



KiX ^ 0) - --Im 



E-iF 



■lEX 



K{X ^ 0) 



8F 



(25) 



-2FX 



where the expression of the positive real numbers E and 
F is given in Eq. (jA6p . 

The transition from one regime to the other takes place 
when two roots of Eq. ([21]) (namely qi and 52) collide on 
the imaginary axis, that is when the discriminant of this 
equation vanishes. This condition yields the expression 
of Afcrit: 



M^. - 1 - -7)2/3 



1 + r;2 + 1 



77^ 



1 



1/3 



1/3- 



(26) 



When 77 — > 0, i.e., in the absence of damping, one recovers 
the usual Landau threshold for emission of Cherenkov ra- 
diation in a weakly interacting Bose gas: Merit = 1 (in 
dimensional units: Vent — Cg). In this case, the perturba- 
tive treatment states that the flow is superfluid for veloc- 
ities below Merit and dissipative above (see Refs. 30 and 
I31I and the computation of the drag in Sec. IIIID|) . This 
is identical to Landau's criterion since both approaches 
give the same value of velocity for the onset of dissipation 
and have the same physical content: excitation of small 
non- localized perturbations is allowed only above Merit- 
In presence of dissipation 7/ 7^ 0, and Eq. (|26|) shows 
that Merit is a decreasing function of 77 (cf. Fig. [2]) . For 
M < Merit (subcritical velocities) there is no Cherenkov 
radiation but, as shown by the explicit computation of 
the drag force below, contrarily to the 7; = case, the 
dissipative effects associated to the finite lifetime of po- 
laritons induce a finite drag force on the obstacle and the 
flow is not superfluid. For M > Merit, Cherenkov radia- 
tion becomes possible but dissipation within the conden- 
sate induces decay of the associated density oscillations. 
The corresponding density patterns are represented in 
each case (M ^ Merit) in the insets of Fig. [2] and the rel- 
evant analytical expressions are given by Eqs. ([24| and 
(ESI). 

The fact that Aferit is modified by damping physically 
explains why perturbation theory is more accurate in 
presence of damping. For a non-damped system, an ob- 
stacle moving at velocity close to Merit = 1 generates Bo- 
goliubov excitations whose typical velocity is also close 
to Cs = 1. As a result, the perturbations accumulate in 
vicinity of the obstacle (since they propagate at the same 
velocity), nonlinear effects cannot be neglected and the 
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FIG. 2. (Color online) Merit = Krit/cs as a function of the di- 
mensionless damping parameter 77, such as given by Eq. (|26|) . 
The dashed lines correspond to the asymptotic expressions 
Merit ==; 1 - |(;7?/2)2/3 (for low Tj) and Merit ==; 2/(3^/377) (for 
large 77). The insets represent typical density profiles in pres- 
ence of a repulsive 5-peak impurity for M < Adait (lower left 
inset) and M > Merit (upper right inset). 



perturbative approach fails^. In presence of damping 
the critical velocity Merit for radiating Cherenkov waves 
differs from the velocity of propagation of small ampli- 
tude perturbation and, moreover, the damping prevents 
large increases of the density. As a result there is no pile 
up of fluctuations in vicinity of the obstacle, nonlinear 
effects may be neglected and the perturbative treatment 
is more likely to be valid. 

This intuitive explanation of the increased accuracy 
of perturbation theory in presence of damping is sus- 
tained by the mathematical reasoning we present now. 
In absence of damping the amplitude of the relative den- 
sity perturbation are of typical mag nitude k/\M'^ -11^/"^ , 
i.e., perturbation theory indeed seriously fails when the 
velocity of the obstacle is close to the speed of sound^^ 
because the expression for 6p diverges. This problem is 
partially cured in presence of damping: for a potential 
of the form ([T9l) a possible estimate of the amplitude of 
\6p{x,t)\ is its value >f |X(0)| at the position of the ob- 
stacle. A study of the dependence of this quantity on the 
velocity and of the damping (i.e., on the dimensionless 
parameters M and 77) shows that, for a flxed value of 77, it 
typically reaches its largest value when M — Merit- The 
value of the quantity >f|ii'(0)| at M = Merit is thus the 
small parameter e of the perturbation expansion, in the 
sense that if this quantity is small for given >c and 7;, the 
perturbation theory is valid for all velocities. This condi- 
tion reads [see formula (IA10|) ] e = - M^^.jJ^/^ < 1. 
Hence e is the small parameter of the perturbation the- 
ory in presence of damping. It never diverges for flnite 
77 and this shows that perturbation theory is more sound 
with than without damping. We see that e effectively de- 
creases in presence of damping because Merit differs from 
1, as advocated in the intuitive discussion of the previ- 
ous paragraph. For small 77, Eq. ([?/)) yields e oc ^irry"^/'^ 



6 



whereas for large rj one finds e oc >r. One can thus equiv- 
alently define the small parameter of the theory as 



e = >f X max{l, 77 ^^^} 



(27) 



and indeed a numerical check shows that, at fixed 77, e 
is a good estimate of the maximum value of for 
X G R and M e M+. 

We stress that the condition e ^ 1 is a criterion of ap- 
plicability of perturbation theory for all M at fixed 77 and 
>f. It is a strong requirement: for given rj and >c failing 
to fulfill the condition e ^ 1, there are still some veloci- 
ties for which perturbation theory holds. For instance in 
the supersonic regime, when T]M{M'^ - l)-^/^ < 1, the 
condition of applicability of perturbation theory relies 
of the smallness of the upstream oscillations and reads 

V(M2 - 1)1/2 ^ ^ 



C. Generic flow pattern for a weak obstacle 

For an obstacle of the generic form (|14p the position 
of the poles of the response function and the critical ve- 
locity ()26p play the same crucial role as for a 5-impurity. 
Eq. ((T5|) yields the following explicit expression for the 
density oscillations: 



X 



Sp{X) = i / dy Res(g3) /cxt(y) e 



/>oo 

i / d?/ V Res(qf)/oxt(y)e' 



qdx-v) 



(28) 



fe{i,2} 



where we recall that Res(g£) is the residue of x(9, —Mq) 
g,t qe {£ = 1, 2 OT 3) [see Eq. Formula ^ is valid 

both below and above A/crit- When 77 = it reduces to 
the one already obtained in Ref.[2^in absence of damping 
[Eq. (45) of this reference] . 

It is interesting to obtain from (pS)) the generic form 
of the long-distance wake which exists ahead of the ob- 
stacle when M > Afcrit- When X is negative and much 
larger than the range of the obstacle potential /ext, the 
first term in ([28)) can be neglected. If, furthermore, /ext 
decreases rapidly enough at —00 so that /oxt(9i,2) ex- 
ists (typically when fcxtix) decreases more rapidly than 
exp[— Im(gi^2) x]), one can approximate the second inte- 
gral by a compact expression yielding 



Sp{X) 



X- 



2Im 



Res(gi)/oxt(gi)e 



iqiX 



(29) 



We recall that Eq. ([^^ is an approximation of formula 
(|28|) valid for M > Afcrit- It is of course exact for 
all X ^ in the case of a (5-impurity. It describes 
Cherenkov oscillations which are damped by a factor 
exp[— Im(gi) x], in complete agreement with the results 
obtained in Ref. [s^ both numerically and also by means 
of Whitham averaging method [Eq. (42) of this reference]. 

Note that for large velocities {M ^ Afcrit) the imag- 
inary parts of qi and q2 tend to zero (cf. Appendix 
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FIG. 3. 5p{X = X + Mt) for a (5-impurity potential (left 
panels) and a Gaussian potential (|30p of width a = 0.5 (right 
panels) as given by perturbation theory [Eq. p5[l ]. The plots 
are drawn for a system in which 77 = 0.5 and in this case 
Merit = 0.5. The two upper panels correspond to a velocity 
below Merit (M = 0.4) and the two lower ones to a velocity 
above Merit (M = 1.75). In the lower right panel the dashed 
gray line correspond to the approximation (|29|) . 



and the wake thus extends far ahead from the ob- 
stacle: the effective damping of the Cherenkov radiation 
tends to zero. However, in this limit, gets very large 
(cf. Appendix E| and for a generic potential |/(qi)| be- 
comes very small: the amplitude of the wake decreases 
uniformly at large velocity, not because of damping, but 
because the large kinetic energy of the flow with respect 
to the obstacle allows to treat this obstacle as a small per- 
turbation. The same effect had been predicted for BEC 
of ultracold vapors in Ref. |2a and has been observed ex- 
perimentally in Refs. ancl|^- 

For being specific, we compare in Fig. [3] the density 
modulations obtained within perturbation theory for a 
(5-impurity obstacle (1191) with the ones corresponding to 
a Gaussian potential of finite width a: 



: exp 



[x + Mtf 



(30) 



When (T — > this potential tends to the (5-impurity po- 
tential (IT^ . As just explained, when M > Afcrit the 
damping of the oscillatory wake in front of the obsta- 
cle is more effective in the Gaussian case than for the 
(5-impurity and is very well described by the asymptotic 
form (P^l as shown in the lower right panel of Fig. [31 



D. Drag force 

In order to discuss the precise influence of the flnite 
lifetime of the polaritons on the possible superfluidity of 
the flow, it is interesting to compute the drag force Fd 
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experienced by the obstacle. Fd is defined as^ 

f dx\,p{x,t)\^d,U,^tix,t). (31) 



A natural way to compute Fd is to insert the perturbative 
expression (IT5|) for 6p in Eq. (PT|) (see, e.g., Ref. Isih . 
Another convenient way is to use the stress tensor T(x, t) 
in a manner similar to what has been done in Ref. l30l . 
The stress tensor is defined as 

T{x,t) ^-Imiri't) + ilV'.r.P - - pC/oxt. (32) 

It verifies the "non-conservation" equation 

^t+T, +p(C/ext)x -2?7(l-p) J, (33) 

where in dimensionless units the momentum current- 
density J coincides with the particle current-density: 
J{x,t) = j{x,t). In presence of damping, in station- 
ary regime, integrating this expression over position, one 
gets 



Fd = 2i^ [ dx (1 - p) J. 

Jr 



(34) 



Within the perturbative approach one can show that 
J{X ^ X + Mt) = -~M5p(X) - 2rjJ^^dy5p{y), and 
using the result ([T5| this yields, for an obstacle of type 



Fd^2r]M [ dx[Sp{x)] 
Jm 

dq 



(35) 



We emphasize that ((3T|) is generally valid, that ([34)1 is 
only valid for a stationary regime in presence of damping 
for an obstacle moving at constant velocity, and that (|35|) 
is the perturbative evaluation of 

For concreteness we now give the explicit expression of 
the perturbative drag (|35p in the case where the potential 
is a Dirac peak of the form ([T^ . One gets 

Fd^ — — 2^ agn{lmqi)qeRes{qe). (36) 

££{1,2,3} 

Substitution of the explicit expressions for the poles 
yields 



Fd 



x2^M(l-M2)-3/2 



cos f (cos f 4- ^ sin |) (cos f + \/3 sin f ) 



3 V— 3 ' V3 3 ' 

for M < Merit and 



Fd = 



F{E^ + 9F2) 



(37) 



(38) 



for M > Merit [in the above expressions 9, E and F are 
given by Eqs. (jAll) and (jA6p ]. The behavior of Fd as a 




FIG. 4. (Color online) Fd/^'^ as a function of M = V/cs 
for different values of the dimensionless damping parameter 
■q. The curves are drawn for the (5-impurity potential (|19|) : 
?7ext(X) = k5{X). 



function of Af is displayed in Fig. |4] for several values of 
rj. For each the critical velocity Aferit is reached exactly 
when the drag is Fd — 2x^/9. The corresponding points 
are shown as white dots in the figure. One can also show 
that for all rj one has Fd = 2yi^ jZ when AI = 1. 
From formulas (ED, (EH) and (EHl one finds 



Fd 



y? X 



r\M when M 
2 when M 



0, 



(39) 



in agreement with the main features of Fig. U) It is in- 
teresting to notice that the drag force is proportional to 
77M when M — >■ (a similar behavior has already been 
observed in Refs. [13 and [27I ). This means that at low 
velocity the obstacle experiences a force which can be 
identified to a viscous drag of Stokes type. When M 
increases and reaches the value M = Aferit, a wake be- 
gins to be emitted ahead of the obstacle. It consists of 
(damped) Cherenkov radiations and one could say, pur- 
suing the analogy with fluid mechanics, that this marks 
the onset of wave resistance. One can push the analogy 
one step further and compare the present results with 
the ones obtained in experimental studies of the drag 
force exerted on objects moving at the surface of several 
viscous fluids. In such experiments it is typically ob- 
served, as in Fig. 21 that the transition to the wave drag 
is continuous^^., but also that Fd considered as a function 
of Y has a quasi-discontinuous behavior for decreasing 
viscosity'^^. An exactly discontinuous behavior is typical 
for the perturbative drag in superfluids'^'^ and is also ex- 
pected on the basis of Raphael-de Gennes theory of wave 
resistance in the context of capillary-gravity waves at the 
surface of inviscid fluids^. This discontinuity disappears 
for finite viscosity^. Moreover, it is interesting to re- 
mark that from Fig. |4] one might erroneously guess (as is 
sometimes done in the analysis of fluid mechanics exper- 
iments) that the relevant critical velocity for the onset of 
wave drag does not depend on viscosity (i.e., on ry in our 
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case) and that at finite viscosity the behavior of Fd{M) 
is just smoothed around the inviscid value [2>c^Q{M — 1) 
in our case]. From our analytical analysis we know that 
in reality the wave drag sets in at Merit [which is not 
equal to the inviscid value Merit (?? = 0) = 1] ^nd that it 
is not possible, when M ~ Afcrit or 1, to disentangle in 
the expression of Fd a viscous component from a wave 
resistance. This is clear from Fig. 2] where the onset of 
wave drag is shown by thick white dots: at these points 
Fd remains a smooth function of M. 

In Fig. m all curves merge at M = 1, and it is in- 
triguing to remark that the drag for a fixed velocity Af 
larger than unity decreases for increased damping. This 
counter-intuitive effect has already been observed in a 
study of the motion of nitrogen drops fioating at the sur- 
face of a liquid batb^. It is explained by the fact that 
viscous effects reduce the range of the wake and accord- 
ingly diminish the wave resistance which is the dominant 
source of drag when M > 1™. 

At large velocity all curves in Fig. 2] tend to the same 
constant value, which is the result for the drag force in 
absence of damping. The fact that the large velocity drag 
does not depend on M is an artifact of the 5-impurity 
potential, as demonstrated by the results obtained in the 
more standard case where the obstacle is described by a 
Gaussian potential of the form ([50)) . In this case formulas 
(|5T|) or ([55]) lead to the expression 

2 

Fd = -Y E 'Z^Res(g,)e-^ 9./2 



sgn(Im qi) + erf 



\V2j 



The corresponding curves are shown in Fig. [S] The 
counter-intuitive 77-dependence already observed in the 
case of a <5-impurity potential is here even more strik- 
ing: the maximum drag is larger at small rj (compare the 
curves obtained for 77 = 0.2 and 77 ~ 0.6). 

In order to better understand the large- velocity behav- 
ior of the perturbative estimate of the drag force we now 
derive an explicit asymptotic expansion valid for any po- 
tential of the form ([Ti)) moving at large velocity. From 
expressions (fT5)) and (|3T|) one gets 



Fd 



dq 
2tt 



qxiq,-Mq)\f,^t{q)\' 
dq 



= -i / dx — qxl?, -Mq)/cxto/oxt(2^)e "^'^ . 



(41) 



In Eq. (|4T|l /cxt ° /cxt is the convolution of /oxt with itself. 
The integral over q in this formula can be evaluated by 
the method of residues. For positive (negative) x the 
contour has to be closed from below (above). Considering 
that when M > Merit the poles qi and q2 which lie in the 
lower half of the complex g-plane verify q2 = —ql and 




FIG. 5. (Color online) Fa/^c^ as a function of M = V/cg for 
different values of the dimensionless damping parameter 77. 
The solid curves are drawn for a Gaussian-impurity potential 
of width a = 0.5. The black dashed line is the corresponding 
asymptotic result (|45[) . The gray dashed line is the result for 
a (5-impurity potential, shown for comparison. 



Res((72) = ~[Res((7i)]*, one gets 
Fd = -2Re 



giRes(gi) / dx /cxt o /oxt (a;) e "^^"^ 

JQ 

qsResiqs) [ da; /oxto/oxt(a:) e"''^"'. 



(42) 



At large velocity one obtains, from Eqs. ([231) ^-^d ([A8| . 

r]M \ 

(43) 



giRes(gi) = -2 + O 
(?3 Resiqs) = O 



(M2 - 1)3/2 



(M2 - 1)3/2 



From this, and using the fact that /oxt ° /cxt is an even 
function of x, one can cast the leading-order contribution 
to Fd in Eq. ((32) under the form 



d^e-iRe(?i).gim(,Obl/^^^o/e,t(a:) 



2 / ^ 



-2Imgi 



R 27r [Re(gi) - q]^+lmqi 



I /oxt (g) I' 



(44) 



The last expression in Eq. ([H[) is obtained using Parse- 
val-Plancherel theorem. At large velocity the imaginary 
part of qi is of order rjM{M'^ — 1)^^, whereas its real 
part is Re(ji ~ = 2(M2 - l)i/2 [cf. Eq. ([M])]: the 
Lorentzian in (j44p is thus a good approximation of the 
Dirac distribution 5{q — qAi)- This directly yields the 
following large velocity result: 



Frf = 2|/ext(9M)|' 



1 + 



■qM 



(M2 - 1)3/2 



(45) 



This means that the typical drag depends on velocity 
(through q^) and tends to zero at large velocity^ con- 
trarily to what occurs for the (^-impurity obstacle. It is 
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interesting to notice that the resuh (1451) does not de- 
pend on 77 at leading order, i.e., that the large- velocity 
drag corresponds to pure wave-resistance. Besides, as al- 
ready remarked in Sec. IIII C[ the obstacle can always be 
treated as a perturbation at large velocity and the asso- 
ciated drag force decreases (the large velocity limit was 
accordingly denoted as "quasi-ideal" in Ref. |42I ) . 



IV. NONLINEAR THEORY FOR A NARROW 
OBSTACLE 

In this section we present results valid for strong ob- 
stacle potentials, in regimes where the perturbative ap- 
proach of the previous section typically fails. In the limit 
of small damping (rj <^ 1) one can expect that other ap- 
proximations are valid. For example, in the case of an 
obstacle represented by a strong (5-potential, one can as- 
sume that the condensate is strongly disturbed at the 
location of the obstacle, so that the difference 1 — p{0) 
is not small; however, the derivative of the distribution 
p{x) downstream the obstacle (for x > 0) is controlled 
by r] and can be considered as small in the case of small 
damping. Hence we can develop for this region a so- 
called hydraulic approximation by neglecting higher or- 
der dispersive effects in our equations (see, e.g., Ref. HI). 
On the other hand, upstream the obstacle (in the region 
X < 0) a supercritical flow generates a stationary oscil- 
latory structure whose oscillation's amplitudes are not 
small, contrarily to what was assumed in the previous 
section. However, in the case of small 77 this oscillatory 
structure can be represented as a slowly modulated non- 
linear wave and, hence, the Whitham modulation the- 
ory can be applied to its description. In this section 
we shall use these two approximate methods (hydraulic 
approximation and Whitham averaging technique) and 
compare their results with the exact numerical solution 
of the problem. 

In all this section we restrict ourselves to the stationary 
version of Eq. ([3]) in presence of a 5-impurity. We find it 
more convenient to work in a reference frame where the 
obstacle is at rest while the condensate moves from left 
to right with an asymptotic velocity and density respec- 
tively equal to M and 1 at both infinities. The equation 
to be solved is the following: 



V 2 



Contrarily to the case of a weak obstacle, where one can 
show that a stationary solution always exists within per- 
turbation theory (see Sec. IIII A| . it is not a priori evident 
that Eq. (j46|) admits a solution. Hence, the assumption 
of existence of a stationary nonlinear regime has to be 
validated by exhibiting the corresponding solution and 
demonstration of its stability. If such a solution cannot 
be found, this means that only time-dependent flows ex- 
ist for the chosen values of 77, >€ and M, which are the 
three parameters characterizing the flow. 



By means of the substitution 

Tpix) = \J p(x) exp i / Ax' u(x') 



(47) 



the Gross-Pitaevskii equation (|46)) can be cast — outside 
the range of action of the obstacle potential — into a hy- 
drodynamic form for the rescaled density p{x) and flow 
velocity uix): 



ipu)x = 2?7/7(l - p), 
^+P 



Px 



Pxx 

1i 



(48) 



We shall use these hydrodynamic notations in this sec- 
tion. 



A. Hydraulic approximation in tlie downstream 
region of a supersonic flow 

In the hydraulic approximation the derivatives are sup- 
posed to be small; hence we can neglect the two last 
terms in the left-hand side of the second of Eqs. PSI) to 
get u2/2 + p = M"^ jl + 1. Then u{x) can be expressed in 
terms of p{x) and substituted into the first of Eqs. 
to give 



P^JM-^ + 2(1 - p<\ =2r^p[\~p) 



(49) 



The solution of this equation, with the boundary condi- 
tion 

p(0) = p, (50) 
can be easily expressed in terms of elementary functions: 




X In ■ 



(1 


-p) 




VI- p + M^/AP H 


h2(l 


-P) 


(1 


-p) 




VI- p + My/M^^ 


h2(l 


-P) 



VM2 



+ = -\'ilj^^^{>i6{x)+ p)i) + ir]{l- p)ij. (46) xln- 



p 


M^ + 2-p+ ^{A'P - 


|-2)(M2- 


h2(l 


-P)) 


p 


AP + 2- p + y/{AP - 


|-2)(M2- 


h2(l 


-P)) 



(51) 

This formula implicitly defines the dependence of the 
density p on x. 

In the supersonic case, in the far downstream region, 
one has 1 — p{x) ^ 1 and one can linearize Eq. (H^ with 
respect to 5p = p — 1. This yields^ 



|(5/9(x)| cx exp 



2r)M 

' AP-1 



(52) 
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The perturbation theory used in the previous section pre- 
dicts the same behavior when r]M{M^ — 1)"'^/^ <^ 1 [Sp is 
found to be proportional to exp{iq3x), where is given 
by (jA8[) ]. However, the range of vahdity of Eq. (|52p is dif- 
ferent: the condition of smahness of the derivative yields 
the following condition of applicability of the hydraulic 
approximation: 



M2 - 1 



< 1. 



(53) 



As a consequence of these different regimes of validity 
one can make the following remark: if 1 — p ^ 1, the 
linearization of Eq. (j49p can be extended down to x = 
0, yielding p(x ^ 0) ~ 1 - (1 - p) exp[-2r] Mx/{M '^ - 
1)]. As we shall see in the numerical section HV C[ this 
approximation has a larger range of validity than the 
pure perturbation approach of Sec. lIIIl This larger range 
of validity of the linearized version of ([5T|) is a result of 
a drawback of the hydraulic approximation: the value 
of p = p{0) is not predicted by this method and has to 
be specified before comparison with numerical results. 
However, we will see in Sec. llVCI that once this is done, 
Eq. ([5T|) gives an excellent account of the downstream 
wave-pattern with slow gradients in a supersonic flow^. 



B. Whitham approximation in the upstream region 
of a supersonic flow 

Upstream the obstacle (when x < 0) supercritical 
flows typically generate a dispersive shock wave which 
is the nonlinear version of the oscillatory wake observed 
in Sec. IHII Now the amplitude of this wave cannot be 
considered as small, but for small rj its parameters are 
poorly modified over one wavelength. Therefore we can 
describe such a flow within Whitham modulation theory 
which is a nonlinear adiabatic approach^. 

The nonlinear progressive periodic wave solution can 
be written in the form (see, e.g., Refs. and 45) 



p{x, t) = i(Ai - A2 - Ag + Xif + (Ai - A2)(A3 - A4) 



X sn' 



(V(Ai - A3)(A2 - A4) (x - t), m) 



and 



u{x, t) =V^ + 



J 



where sn is the sine elliptic Jacobi function, 

(Ai-A2)(A3-A4) 



(Ai- A3)(A2- A4)' 



(55) 



(56) 



The parameters Ai ^ A2 ^ A3 ^ A4 are called the Rie- 
mann invariants of the system. In the case of strictly 
periodic solutions they are constant and they determine 
characteristics of the wave such as the phase velocity V^p 
[Eq. (j56p ]. the current j evaluated in the frame where the 
wave is standing [Eq. ((57)) ]. the amplitude of the oscilla- 
tions 



a = (Ai - A2)(A3 - A4), 
and their wavelength 



L = 



2K(m) 



V(Ai-A3)(A2-A4)' 



(58) 



(59) 



K(m) being the complete elliptic integral of the first kind. 
In the modulated dispersive shock wave occurring in the 
upstream region, the A's become functions of x which 
vary weakly over one wavelength. We consider here the 
stationary solution and hence these parameters do not 
depend on time t and the phase velocity is equal to 
zero: 



0. 



(60) 



However, in the upstream region, the A's are functions 
of position and their x-dependence is determined by the 
Whitham equations (see Appendix |B]) 



dA^ 
dx 



2 GiAi + G2 



{1,2,3,4}, (61) 



where 



Gi 
G2 



-ri I Av , 



V) 



' av ■ 



(62) 



TZ{v) and vi, V2, 1^3 being defined by Eqs. (jB14|) and 
(IB15[) . According to Eq. ([TO)) , the system (pT|) admits 
(54) the first integral = 0. We shall now show that 

it admits another integral and can thus be reduced to a 
set of two (coupled) differential equations. To this end, 
we shall use the Jacobi identities 



4 

E 



A* 



~[ Wm^i'yXi ^ Am) 



for ^ fc ^ 2, 



and 



4 

E 



(63) 



4 Wm^ik^i ^ Am) 



and 



j = |(-Ai-A2 + A3 + A4) 
X (-Ai + A2 - A3 + A4)(Ai - A2 - A3 + A4). 



(57) 



to obtain #i = #a and 

da: da; 



ds3 
dx 



2Gi 



ds4 
dx 



2^2 



(64) 
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where the s's are symmetric functions of the A's: 

i<3 



S2 — ^i^Ji 



S3 



i<j<k 



K^j^ki S4 — A1A2A3A4 



(65) 



Here si and S2 are the integrals of our system. The value 
of Si is already known from Eq. (160^ : si = 0. In order 
to determine the value of S2 we calculate the asymptotic 
values of the Riemann invariants at x — ^ —00, where the 
flow is stationary with p = po = 1 and u = uq = M > 0. 
The amplitude of the oscillations vanishes here; hence we 
find from (|58p that Ai = A2 (another possible choice is 
•^3 = -^4; it corresponds to a flow with M < 0). Then, 
from Eq. (jBlSp we have the equation 



(66) 



lim p{x) ^ pq^I 

X—^ — OO 

= vi = V2 = i(A3 - Xif, 
as well as the expression for the current density 
lim j{x) = pqUq = M 

= |(A3-A4)2(-2Ai+A3 + A4), 

from which we get another equation: 

M= i(-2Ai + A3 + A4). (68) 

With account of Eq. (jMl) (that is 2Ai + A3 + A4 = 0) we 
find, at X —00, 



(67) 



Ai = A2 
Hence, 



M , M , M 

y, A3 = y-1, A4 = y + 1. (69) 



/f2 



(70) 



As a result, we can define the functions Ai — Ai (53,54) 
as being the roots of the equation 



A* 



1 A^ - 53 A + 54 = 0, 



(71) 



ordered according to Ai ^ A2 ^ A3 ^ A4. Substitution 
of these functions into (|B15P and of the results into ((64)) 
yields the system of two differential equations for 53 and 

54, 



d53 _ 2Gi(53, 54) d54 _ 26*2(53,54) 



da; 



L(53,54) 



da; 



L(53, 54) 



(72) 



We now have to find the initial conditions for this system, 
that is to determine the values of 53 and 54 at a; = 0. To 
this end, we take into account that Whitham theory im- 
plies that the parameters of the wave weakly change over 
a distance of about one wavelength. Therefore we can 



assume that, to the left of the obstacle and close enough 
to it, the wave can be approximated by the cnoidal wave 
solution (|54|) . (|55|) and to the right of the obstacle it is 
given by a hydraulic approximation parameterized by the 
the value p of the density at the location of the (5-obstacle. 

It is known (see, e.g., Ref. E3) that a non- modulated 
cnoidal wave solution p(x) satisfies the equation 



p, = 2^n{pj, 

where the coefficients of the polynomial 
7^(^y) = ~ - i'2){v - vz) 
= v'^ + 2s2v'^ + [si ~ 454)1/ 



(73) 



(74) 



are expressed in terms of the symmetric functions 52, 
S3 and 54. Then the solution ([M]), ([55]) (with = 0) 
can be expressed in terms of the zeros ui, 1^2, 1^3 of this 
polynomial as follows 



u{x) 
where 



p[x) = vi + [u2 — 1^1 ) sn^ (Vt^3 — ^1 X, m) , 
S3 



p{x)' 



m 



V2 - v\ 
vz - vx 



and L 



2K(m) 
- V\ 



(75) 



(76) 



In the stationary modulated situation we consider that 
v\,V2.,V'i, sz,rn and L do not depend on time in Eqs. (|75|) 
and (|76p . but they all depend on x. 

It follows from Eq. (|17p that the current of polaritons is 
preserved in transition through the (5-potential; j (O^) = 
j^O"*"). Then the second of Eqs. ([75]) yields (under the 
assumption that the hydraulic approximation is valid for 
a; ^ because <C 1) the value of 53(0): 



S3 



(0) = ^.(O)p(O) = -p^M-^^KX-p). (77) 



For calculating the value of 54(0) we use the matching 
condition at x = 0: 



P.(0+)-p.(0') =4xp(0). 



(78) 



Pursuing the use of the downstream hydraulic approxi- 
mation already used in Eq. (j77l) we write p(0) = p and, 
from Eq. (gS]), 



P=.(0+) = 



277p(l-p)v/M2 + 2(l-p) 
A/2 + 2 - 3p 



(79) 



In the same spirit of a small-77 approximation we have 
from Eq. ^ Px[^^) = -2^n[p), so that Eq. ^ reads 



nip) 



2>cp ^ 



Vp{l~p)^AP + 2(l-p) 



M2. 



3p 



This yields 

(1/11/2 + iy2'^3 + i's.t^i)x=o = si - 454(0) 
_ [A^p~pxiO+)r 



Ap 



{2A'r 



(80) 



(81) 
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and then 



S4(0) = - 



(M2 + 2)2 



2(l'r + 2)p + 3p' 



2{AP + 2 - 3p) 



(82) 



We note here that for smaU values of rj, the analytical 
expression (|82p can be simplified by replacing Eq. ((79)) 
by the simple approximation Px(O^) = 0. This amounts 
to also replace 77 by in the expressions (|80|) and (|82|) . 
This simple scheme is accurate when r/ < 0.5. 

Equations ((77)) and ((82p give the initial conditions for 
the system 



dS3 

dx 

ds4 
da; 



277 
"IT 

L 



dz^ — , 

''1(S3,S4) V7^(^^) ' 

l^2{s3,Si) 1 — 1/ 

dz/ , 



(83) 



where ^4(53,54) (i = 1,2,3) are determined as being the 
roots of the equation TZ{p, 83,34) = 0, where 



TZ{iy, S3, S4) = ly^ 



(Af2 + 2)^ 
2 



(84) 



In Eqs. (l83l) L is also expressed in terms of the i^'s [see 

Eqs. dzg)]. 

In the present application of Whitham modulation the- 
ory it is important to notice that for fixed values of >c, 
77 and Af , the solution of Whitham equations depends 
on a single parameter p which is also a function of the 
same set of physical parameters {k, j], M) prescribed by 
the external potential and the boundary conditions of 
the Gross-Pitaevskii equation. Hence, the parameter p 
can be found from the condition that the solution of 
Whitham equations satisfies the correct boundary con- 
dition at X — >■ —00, namely that the envelopes of the 
density oscillations tend to the asymptotic value of the 
density: 



vi{x), 1/2 (x) 1 



(85) 



Some values of p calculated in this way are listed in the 
second row of Table U for r] — 0.05, M = 3 and several 
values of x. We compare them with the values of p ob- 
tained by exact numerical solution of Eq. ((^S)) . As we 
see, the agreement is very good. 

The dependence of p on 7; is displayed in Fig. |6] (left 
panel) for several values of x. This plot suggests that 
p — p{Q) does not tend to unity in the limit 77 — ;> 0. This 
means that, in this limit, the flow pattern does not re- 
duce to the exact solution found in Ref. in the case 
77 = 0, since for this solution p(Q) = 1. Rather, however 
small is 77, 1 — p remains finite, the wave structure occu- 
pies a portion of space proportional to r/~^ and decays 
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0.5056 


0.4011 


p (numerics) 


0.9352 


0.7916 


0.6377 


0.5055 


0.4013 



TABLE I. Values of p for different values of >c in the case 
M = 3 and 77 = 0.05. The row p (Whitham) corresponds to 
the value of p = p(0) found by solving Whitham equations 
(|83|) and imposing the condition (|85|) (see the text). The row 
p (numerics) corresponds to the value of p(0) found via a 
numerical resolution of Eq. (|46|l . 
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FIG. 6. p = p(0) as a function of rj for M = 3 and different 
values of h (left panel) and as a function of x for ij — 0.05 
and several values of M (right panel). 



towards an undisturbed fiow (p = 1) at 3> 77^^. The 
dependence of p on >f for several values of M is shown in 
the right panel of Fig. (6] 

A striking feature of the plot in the left panel of Fig.jHjis 
the extremely weak 77-dependence of p. This important 
property of the theory can be explained by the simple 
fact that the space coordinate x and the parameter 77 en- 
ter into both the hydraulic approximation and Whitham 
equations only through the combination 77 x [see Eqs. ((5ip 
and ()83p ]. As a result, 77 can be rescaled out of the ex- 
act relation ((T5)) after averaging over fast oscillations in 
the dispersive shock region a; < 0, so that we arrive to 
an equation which depends on 77 only through the small 
value of p:r(0+) [see Eqs. ((TD) and ^]. If we neglect this 
term, then the resulting equation yields p as a function 
of M and x only. 

When p is found, all the parameters of the dispersive 
shock wave are determined, the functions 7^i(x), i'2{x), 
1/3(2:) can be computed by solving Eq. ()83p . and their 
substitution into Eq. ((75)) yields the oscillatory structure 
upstream the obstacle. The same value of p determines 
the hydraulic solution downstream the obstacle. Thus, 
we reach a complete description of the nonlinear wave 
generated by a supercritical flow past a (5-obstacle. 

The accuracy of the theory is illustrated by Fig. [7| As 
we see, the agreement between the results of the com- 
bined Whitham and hydraulic approaches and the nu- 
merical computations is excellent. Note that Whitham 
method is perfectly valid in a regime where the pertur- 
bative theory of Sec. IIIII seriouslv fails (|p(a;) — 1| is not 
small in Fig. [7]). For illustrative reasons we have chosen a 



13 




— - Numerics 

— Whitham theory 

— Hydraulic approximation 



-15 



-10 



-5 



10 



X 



FIG. 7. (Color online) Comparison of the Whitham theory 
with the numerical solution of Eq. (|46p . The plot is drawn in 
the case rj — 1, M = 3 and x = 4. The numerics corresponds 
to the dashed black line. Whitham envelopes are shown by 
thin red solid lines, and the upstream dispersive shock wave 
oscillatory structure obtained by substitution of the solution 
of the Whitham equations (|83p into Eq. (|75p is shown by a 
red solid line (for x ^ 0). The downstream (x ^ 0) hydraulic 
approximation is shown by a green solid line. 



relatively large value of ry {rj = 1): we wanted to work in 
a regime where the overall modulations of the oscillating 
pattern occur over a characteristic length which is not to 
large with respect to the wavelength of the oscillations. 
As we see, even in this unfavorable case the agreement 
with the exact numerical results is very good. 

The solution of the system (|83)) exists, and the up- 
stream pattern can be described as a slowly modulated 
cnoidal wave, as long as its initial conditions can be 
found, that is, as long as the equation Tl{h', 33(0), S4(0)) = 
has three real roots. If 77 is strictly zero, then p — 1, 
and this equation reads 

+ (M^ + 2)z/2 + (1 + 2Af2 + 4>f2)i^ - = 0. (86) 

Two of the roots coalesce and go into the complex plane 
when the discriminant of Eq. (I86p vanishes. This cor- 
responds to a boundary between possible parameters in 
the plane (>f, Af ) determined by the condition 



1 

32 



(87) 



The same boundary was already found in a different an- 
alytical form in Ref. 29 for a non-damped system. In our 
problem (r/ 7^ 0, p 7^ 1), this boundary is changed and 
can be found by numerically determining when the dis- 
criminant of Eq. ()84[) vanishes. However, when 77 7^ 0, as 
we shall see in Sec. lIV Cl new stationary solutions appear 
when K gets so large that the upstream flow is not de- 
scribed by a modulated cnoidal wave, making the deter- 
mination of the domain of validity of Whitham approach 
less crucial than when 77 = 0. 
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FIG. 8. (Color online) Different profiles p{x) for flows past 
a (5-impurity potential of type (I19|l . For all the profiles the 
damping parameter is rj — 0.05. For the upper row M = 3, 
then M — 1.2 for the row below, M = 1 for the following 
one and finally M = 0.5 for the lower row. The value of >c 
is indicated in each plot. In each plot the black solid line 
corresponds to the numerical solution of Eq. (|46|l . the (red 
online) thin line to the perturbative result and the (green on- 
line) dashed line to the result of the hydraulic approximation 
which is only relevant for x ^ (see Sec. IIV A| . 



C. Numerical results 

In this section we present results of the full numerical 
solution of Eq. (^5)) . We used a shooting method, starting 
the numerical integration from large and positive x with 
an initial behavior given by the prediction of perturbation 
theory. Typical results are displayed in Fig. [8] 

The upper plots of this figure are drawn for M = 3 
which is a velocity deep enough in the supersonic regime 
for Whitham theory of Sec. IIVBI to apply over a rather 
large range of values of x. The left plot of the upper row 
corresponds to >f = 0.5. For this value of perturbation 
theory is valid upstream (x < 0) but fails for positive x, 
whereas the hydraulic approximation is quite accurate in 
this region, as shown by the dashed line in this plot. For 
X- = 4 (right plot of the upper row of Fig. the density 
profile shows the same features, but in this case perturba- 
tion theory seriously fails, whereas the downstream wave 
pattern being typical for a damped cnoidal wave is very 
well described by Whitham theory (not shown, because 
undistinguishable from the numerical result). 

The two rows below the upper one correspond to 
M = 1.2 and M — 1. They are interesting because they 
show that, whereas perturbation theory fails in absence of 
damping when Af ~ 1, for 77 7^ it has a regime of valid- 
ity even for velocities M close to unity. This is illustrated 
by the good agreement of the perturbative results with 
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the numerics displayed in the two left plots of the central 
rows (which are both drawn in the case >c — 0.05). It is 
also interesting to remark that for M = 1, no stationary 
solution exists when 77 = 0, whereas here we could find 
such solutions up to >tf = 0.3 (see the right plot of the 
third row). The values M — 1 and x — 0.3 are close to 
the boundary marking the end the existence of station- 
ary solutions when 77 = 0.05. In this case we see that the 
downstream wave pattern shows small scale disturbances 
which were recognized in Ref. [s^ as typically occurring 
near the end of the stationary regime. 

The second upper row of Fig. [5] corresponds to M = 
1.2. In this case, when ij = 0, there is no stationary 
solution for >c > 0.0495 [see Eq. ^ or Ref. HI]. As seen 
on the figure, when 77 = 0.05, one can find stationary 
solutions for much larger values of >c (up to >c ~ 1.2; 
see the corresponding plot). However, the density profile 
found in this case is very different from a damped cnoidal 
wave. It seems to be a stationary version of a type of 
time-dependent profiles studied in Ref. 46 for the case 
?7 = 0: a plateau develops just upstream the obstacle 
which terminates when x — ?> —00 by a dispersive shock 
wave. Here, when 77 7^ 0, the plateau and the shock wave 
are damped because the specific form of the modified 
Gross-Pitaevskii equation ^ favors relaxation towards 
p=l. 

The lower row of Fig. [S] displays results corresponding 
to a subsonic obstacle {M = 0.5). For this value of the 
velocity, there is no stationary solution in the 77 = case 
for X ^ 0.5 9^^'"^^ . As illustrated by the right plot of this 
row (drawn for x = 0.7) in presence of damping, solu- 
tions exist for slightly larger values of x. However we find 
that, when 77 passes from to 0.05, the range of values 
of >f allowing for a stationary solution does not increase 
in the subsonic case as much as it does in the supersonic 
region. This is illustrated by Fig. |9] were we represent 
the domain of existence of stationary flows in the (x, M) 
plane. This domain corresponds to the shaded region in 
the Figure and was numerically determined in the case 
77 — 0.05. The large increase of the stationary domain 
for supersonic flows in presence of damping is due to 
the occurrence, when 77 ^ 0, of a new class of profiles 
with an upstream plateau, as explained above and illus- 
trated in the right plot of the second row from the left 
in Fig. |H] (corresponding to Af = 1.2 and x ^ 1.2). This 
type of profile cannot be stationary with the boundary 
condition p{x — ±c») = 1 in a non-dissipative system. 
Here, the damping term in Eq. ([3|) provides a mecha- 
nism allowing the downstream relaxation from p(0) < 1 
to p{x 00) — 1 and the upstream dispersive shock is 
stabilized by dissipation. 

The inset in Fig. |9] represents the exact domain of sta- 
tionary flows for 77 = 0, as analytically determined in 
Ref. [2^. One can identify three regimes depending on 
the value of the parameters x and M = V/Cs'- (i) sub- 
sonic, stationary and superfluid, (ii) dissipative and time- 
dependent, (iii) dissipative, stationary and supersonic. 
As seen in this inset, regimes (i) and (iii) are always sep- 
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FIG. 9. (Color online) Different regimes of flow past a 5- 
impurity potential of type (|19p in the (>f, M = V/cs) plane. 
The main plot corresponds to the dissipative system (with 
77 = 0.05). The inset is drawn for — 0. In both plots the 
shaded regions correspond to stationary flows, the white ones 
to time-dependent and dissipative flows and the horizontal 
dashed line indicated the transition from a localized wake to 
a regime of Cherenkov emission as predicted by perturbation 
theory. In the inset, the subsonic (M < 1) shaded region 
is superfluid and the supersonic (M > 1) one is dissipative; 
the exact equation of the boundaries between the different 
domains is given in Ref. [2^. In the main plot the dots with 
error bars represent the numerically determined boundary of 
the stationary domain. They are connected by a dashed line 
to guide the eye. The other dashed line in this plot represents 
the rj = result (shown for comparison). 



arated by the time-dependent region (ii). This feature is 
also valid for a thick obstacle'^'* and is in contradiction 
with the (wrong) prediction of perturbation theory for 
77 = 0. Indeed, in the non-dissipative case, perturbation 
theory always fails when V is close to and in this case 
the true flow gets time-dependent. On the contrary, for 
finite ^ 1 we showed in Sec. Illll that the perturbative 
prediction of existence of a stationary flow pattern for all 
velocities is valid until x ^ rj^^^ . This is corroborated by 
the numerical results displayed in Fig. IHlfor 77 — 0.05. In 
this case rj^^^ ~ 0.3 whereas the largest value of x for 
which a stationary flow exists for all M is numerically 
found to be ~ 0.1. Then, we can go one step further, 
and since we showed that the actual small parameter of 
perturbation theory is e — x x max{l, 77^^/'^} we con- 
jecture that the neck of the stationary domain in Fig. m 
extends when 77 increases from 0, until 77 ^ 1, where the 
largest value of x for which a stationary flow exists for all 
M should remain approximatively constant and of order 
of 1. 



V. CONCLUSION 

In the present work we have analyzed the flow of a one- 
dimensional polariton condensate in motion with respect 
to an obstacle in a situation of non-resonant pumping. 
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We solved the problem perturbatively and showed that 
at this level there exists a smooth crossover from a viscous 
flow to a regime where the drag is mainly dominated by 
wave resistance. Perturbation theory predicts that this 
occurs at a velocity Merit independent of the potential 
representing the obstacle. We argued that in the case 
of a (5-impurity [represented by a potential of type (1191) ] 
the perturbative approach is valid for all velocities in the 
regime >f x max{l, ry~^^^} <C 1, where rj is the dimension- 
less damping parameter defined in Eq. 1^. As shown in 
the previous section this implies that stationary profiles 
indeed exist for all velocities ii >c < min{l, yy^/'^}. In this 
case there is a continuous transition from a dissipative 
drag to a regime dominated by the wave resistance. 

However, from Fig. |9] we are led to refine this discus- 
sion of the transition between a regime where the wake 
is localized in vicinity of the obstacle and a regime of 
(damped) Cherenkov radiation: we see on the example 
of the (^-impurity that for a strong enough potential the 
two types of flows are separated by a time-dependent 
regime, as typically observed in BEC atomic vapors. In 
this case one cannot state that the crossover is smooth. 

An important result of our work is the demonstration 
that it is difficult to assess on the superfluidity of a po- 
lariton system just by studying the density perturbation 
past a localized obstacle. In particular, we showed that 
the absence of long-range wake cannot be used as a cri- 
terion for the absence of dissipation. 

The comparison of our results with the ones of Ref. [13 
leads to the conclusion that the gross features of the wave 
pattern discussed in the present work are quite indepen- 
dent of the technique used for setting the fluid into mo- 
tion with respect to the obstacle. However, we use a spe- 
cific model [Eq. ([1])] with non-resonant pumping which is 
more relevant for the experiment presented in Ref. 9. In 
this experiment, a two-dimensional supersonic cloud of 
polaritons colliding with an obstacle was observed to in- 
duce a rather well defined wake, with oscillations having 
an apparently specified wavelength. The same feature 
was observed numerically in Ref. (see also the dis- 
cussion in Ref. IH). Our perturbative results allow to 
understand this phenomenon in a one-dimensional set- 
ting: the pattern of the upstream oscillatory wake in a 
supercritical fiow {V > Vent) is governed by the complex 
wave vectors qi and q2', see Sec. lIIII Also in the nonlinear 
approach (Whitham theory of Sec. lIVBp the wake keep 
a simple shape: perturbation theory fails to properly ac- 
count for the amplitude of the oscillations, but it still 
approximatively describes their wavelength. 

Finally, this work naturally calls for developments. 
One would first like to precisely determine the domain 
of time-dependent nonlinear flows in presence of damp- 
ing. Secondly, one would like to extend the present work 
for taking into account polarization effects, and, thirdly, 
it is natural to apply the perturbative approach to higher 
dimensions. Works in these directions are in progress. 
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Appendix A: Poles of the response function 

X(q, -Mq) 

In this appendix we determine — as a function of M — 
the location in the complex g-plane of the poles of the 
response function evaluated at w = —Mq. Consid- 
ering the expression of x one sees that these poles are 
the three zeros of D{q, —Mq)/q. We denote them as gi, 
92 and q^. They are solutions of Eq. ([2T|) . This equa- 
tion has three imaginary solutions when its discriminant 
A = 256(1 - APf/27 - MrfM"^ is positive. The condi- 
tion A > is equivalent to M < Merit where the expres- 
sion of Merit is given in Eq. (|26)) . In this case, deflning 



/877M\ 
arctan — 

V Va; 



(Al) 



one finds 



;i-M2 fe TT 
'Zi=4:\/^^sm(--- 



q2 = -4i \l sm g 



,3=4M/^^sm - + - 



Alternatively one can write qi = i{—A + B), q2 
and qs = i{A + B) with 



cos(6l/3) 
^sin(0/3) 



= 2a/i - M2 
If A < 0, i.e., if M > Merit, defining 



1/3 



(A2) 



-2\B 



(A3) 



(A4) 



one finds 

qi — D(^^-j exp(— i7r/6) — D(-) exp(i7r/6), 
92 = cxp(i7r/6) + cxp(-i7r/6), 

g3 = i(i^(+)+^(-)). 



(A5) 



Alternatively one can write qi = E — iF, q2 = —E — \F 
and (73 — 2iF with 



V3, 



1 



E=^{D(+)-D(_)), F^-{D(+)+D(_)). (A6) 
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FIG. 10. (Color online) Position of qi, q2 and 53 in the com- 
plex g-plane. The figure is drawn in the case r] — 0.1. The 
arrows indicate the direction of motion of the poles when M 
increases from to 00^. 



One can verify that X^fci Qt — for all values of M, as 
already clear from the form of Eq. (HI]). A similar relation 
holds for the residues of —Mq) whose expressions are 

given in ([231): J2e=i ^^^{qi) = 0. 

The typical M-dependence of the position of the poles 
in the complex plane is illustrated in Fig. [TO] When 
M = one has 9 = 0, q2 ~ and (73 = —qi = 2i. When 
M is increased from zero, qi and 52 get closer on the 
imaginary axis until they collide (when M = Merit) and 
then acquire a finite real part. When M — > 00, 53 — > i 0+ 
and 1 -J — >■ (±)oo — 10+. 

A useful approximation for the expression of the poles 
is obtained when rjM/lAP — V^l"^ <^ 1. In this case one 
obtains, when M < Merit, 
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(zp)2Vl-M2 + 
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(A7) 



and when M > AL- 



(7(1) ~ (±)2Vm2 - 1 - i 
rjM 



T]M 

M2 - 1 ' 



(A8) 
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The above expressions are valid up to corrections of rela- 
tive order ?72m2/|M2 — Ip. It is interesting to notice that 
expansions (IA7p and (|A8p are equally valid at large ve- 
locity and at small damping. Indeed, as discussed at the 
end of Sec. IIII Dl at large velocity the effects of damping 
are negligible. 

From the explicit expressions (jA2l) and (|A5I) of the 's 
it is a simple matter to evaluate the integral (jl6p which 
permits to compute the function K{X). One gets 



K{X ^ 0) 

K{x s; 0) 



iRes(g3)e'«-'^, 

-i [Res(qi) e'91-^ + Res(q2) e'«^-^] 



(A9) 



Formulas (|A9[) are valid for all M, but the explicit ex- 
pressions for the qis depend on M. For instance, when 
M < Afcrit the q^'s are all imaginary and K tends 
rapidly to zero when |X| — >■ 00. On the other hand, 
when M > Merit the exponential decrease of K{X) gets 
weaker (because the imaginary part of the q^'s is smaller) 
and K{X ^ 0) oscillates (because qi and q2 acquire a 
real part). The typical density perturbations associated 
with K [i.e., for a S-peak potential of the form (IT9l) ] are 
sketched in the insets of Fig. [21 Note that the value of the 
qi's does not depend on 77 when M = 0, i.e., within the 
theoretical description corresponding to Eq. ([Ij , the den- 
sity perturbation induced by a motionless obstacle does 
not depend on the damping. 

The expressions (jA9l) are equally valid in absence of 
damping, i.e., when 77 = 0. In this case Afcrit = 1, 92 = 
and (73 = —qi for AI < Merit and for M > Merit, 93 = 
whereas qi and (72 are real and opposite [cf. Eqs. (IA7[) 
and (|A8[) ]: for M > Merit and rj — one observes un- 
damped Cherenkov radiations ahead of the obstacle as 
discussed in the main text. For M > Merit and r? 7^ 
these Cherenkov radiations are damped since in this case 
qi and (72 have a nonzero imaginary part. 

Finally, we need to evaluate the order of magnitude 
of the quantity 1/^(0)1 at AI = Aferit since, as argued 
in the main text (Sec. IIII Bp . this is the small param- 
eter of perturbation theory for a (5-impurity obstacle. 
For M = Merit one gets 9 = 7r/2, qi = q2 = -93/2 = 



-2i^(l - Af2.J/3 [cf. Eqs. (|A2l)] and this yields 



xK{Q) 



when M = Mr 



(AlO) 



From the expression (|26p for Aferit one sees that (1 — 
-^crit)"^^^ - 75(2/??)^^^ when 77 < 1 and tends to unity 
at large M, from which one obtains the estimate (|27p. 



Appendix B: Derivation of perturbed Whitham 
equations 

The general method of derivation of the Whitham 
equations for perturbed integrable equations which in 
their nonperturbed form belong to the Ablowitz-Kaup- 
Newell-Segur scheme was developed in Ref . 113 and it can 
be formulated as follows. Let the evolution equations of 
some field variables Uk have the form 



dt 



dUr. 



dx 
dx 



9a;2 ' 



i9a;2 



(Bl) 



where a small parameter e ^ 1 is introduced which mea- 
sures the dispersion effects. It is supposed that a non- 
perturbed system 



duk 
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dx 
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(B2) 
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can be represented as a compatibility condition of two 
linear equations 



Ax, 

Xt = -^BxX + ISxx, 



(B3) 



where A and B depend on the Wfc's, their space derivatives 
and on the spectral parameter A. It is assumed that 
the system (|B2p has a periodic solution with wavelength 
L (X e and it is parametrized by the constant parameters 
Ai which appear in the finite-gap integration method in 
the following way. The second-order linear equation (jB3p 
has two basis solutions x± and their product g = x+X- 
satisfies a third-order differential equation which can be 
integrated once to give 



-j99xx 



aPiX), 



(B4) 



where a is determined by the sign of the highest order 
term in ^ as a function of A, i.e., A ~ —aX" as A — )■ cx). 
Periodic solutions are distinguished by the condition that 
P(A) is a polynomial in A and then A^ are its zeros. We 
shall confine ourselves to the one-phase periodic solutions 
which physical variables depend on a single variable x — 
Vip t only. 



In a modulated wave the parameters Ai become slow 
functions of x and t whose evolution is described by the 
Whitham equations which in the case (IBip . (IB3p can be 
written in the form 



^A, {B/g) d\, 
dt (l/g) dx 




1/5) n™^.(A.-A. 




dA d^^'-'^Rk 




(B5) 



where ik denotes the highest order of derivative of Uk 
entering in A. The angle brackets denote the averaging 
over one wavelength: 



m = 



dx T . 



(B6) 



The spectral parameter A should be put equal to Ai after 
averaging. 

We shall apply here this scheme to the perturbed non- 
linear Schrodinger (NLS) equation 

i£ i,t + le'^xx - 1^1 V = iG(|^n^, (B7) 

where G'(p) is a real function of the density p — 
Eq. (I3|) pertains to this type [with G(p) — 77(1 — p)]. In 
the case of Eq. (|B7|) we have two field variables "0, i/'*! 
and, correspondingly, two terms of perturbation in (|B1I) : 



For non-perturbed NLS equation the linear system (jB3 
is specified as 



(B8) 



A = -\' -\e\^-.r<^~^^ + -r^, (B9) 
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Substitution of (IB9p into (IB5P shows that, in the expres- 
sion to be averaged [in the right-hand side of (jBSP J. the 
leading term in powers of e is equal to 2{Gpg)/e. The 
averaging can be performed with the use of equations 



known from the theory of periodic solutions of the NLS 
equation (see, e.g., Ref . lish : 
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where P{p,a) = Y{i{^J'a - K) and si = Ai- The quan- 
tity Pa is known as the auxiliary eigenvalue in the finite- 
gap integration method. Hence, we obtain 



9 / \ A - Aia 

g / L dXi 



2 dL 



(B12) 



For calculating (Gpg) we also take into account that pa 
can be expressed as a function of p in the following way 
(see Ref. Hi): 



Ma(p) 



4 



-J +1 



2p 



(B13) 



where 



7^(^y) = (zy-i/i)(jy-jy2)(i/~i/3), ^^^^^ 
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J^i - i(Ai - A2 - A3 + A4)^ ^jjg stationary case, i.e., when d\/dt = and si 

1^2 = |(Ai — A2 + A3 — X^y, (B15) 2V^ — 0, the Whitham equations simphfy to 

l^3 = 3(Al+A2-A3-A4)^ 



and 



e^ = 2Vn. (B16) dA, _ 2 dA, + G2 . 

dx -Ln™^.(A.-A„)' ^^^'^ 



da; 

Then we obtain the Whitham equations for the Riemann 
invariants A^ in the form 

dXi , dXi Vi - si/2 

'u - 



dt dx ]\r,l^^{X^ " Xm) where 

2 n G{v)[{\-s,/A)u~j/2] 
X— / dv , 

(B17) 

with Gi=-/ dz.-^, G2 = -'- Ay 



(B20) 

= "T+(t^) ' «e {1,2,3,4}. (B18) For G'(p) = 77(1 - p) we arrive at Eqs. dH]) and dM])- 
L oXi 



si f2dL 
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